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Abstract 


We  study  a diffusive  energy-balance  climate  model, 
governed  by  a nonlinear  parabolic  partial  differential 
equation.  Three  positive  steady-state  solutions  of  this 
equation  are  found;  they  correspond  to  three  possible 
climates  of  our  planet:  an  interglacial  (nearly  identical 

to  the  present  climate),  a glacial,  and  a completely  ice- 
covered  earth.  We  consider  also  models  similar  to  the  main 
one  studied,  and  determine  the  number  of  their  steady  states. 
All  the  models  have  albedo  continuously  varying  with  latitude 
and  temperature,  and  entirely  diffusive  horizontal  heat 
transfer.  The  diffusion  is  taken  to  be  nonlinear  as  well  as 
linear . 

Wo  investigate  the  stability  under  small  perturbations 
of  the  main  model's  climates.  A stability  criterion  is 
derived,  and  its  application  shows  that  the  "present  climate" 
and  the  "deep  freeze"  are  stable,  whereas  the  model's  glacial 
is  unstable.  A variational  principle  is  introduced  to  confirm 
the  results  of  this  stability  analysis. 

We  examine  the  dependence  of  the  number  of  steady  states 
and  of  their  stability  on  the  average  solar  radiation.  The 
main  result  is  that  for  a sufficient  decrease  in  solar  radia- 
tion (about  2 percent)  the  glacial  and  interglacial  solutions 
disappear,  leaving  the  ice-covered  earth  as  the  only  possible 


climate . 


1.  Introduction 


The  concept  of  a climate  is  one  of  those  abstractions 
which  appears  to  be  self-evident  to  the  layman,  but  is  by 
no  means  well  defined  scientifically.  The  intuitive  idea 
of  a climate  has  two  aspects: 

(a)  the  most  important  features  of  atmospheric  phenomena; 

(b)  the  average  behavior  of  these  phenomena  over  a suitably 
long  time  interval  and  over  sufficiently  large  areas. 

The  difficulties  start  when  one  tries  to  give  a precise 

meaning  to  the  key  words  "most  important",  "suitably  long" 
and  "suitably  large".  We  start  with  "suitably  long"; 
clearly,  a year  is  an  absolute  lower  bound  for  a reasonable 
averaging  time  interval,  since  daily  and  seasonal  variations 
should  be  excluded.  To  decide  over  how  much  longer  than  a 
year  the  averaging  should  be  performed,  one  has  to  look  at 
the  record.  There  are  three  kinds  of  records:  instrumental, 

the  length  of  which  is  of  the  order  of  hundreds  of  years, 
historical,  of  the  order  of  thousands  of  years,  and  geological, 
of  the  order  af  hundreds  of  thousands  of  years.  These 
records  show  that  features  of  the  atmosphere  change  on  all 
the  time  scales  represented  in  them  (e . g. , Robinson,  1971). 

Thus  it  would  appear  at  first  that  it  is  not  possible  to 
distinguish  between  "fast"  variations  in  yearly  averages, 
which  should  be  averaged  out  when  defining  a climate,  and 
"slow"  variations,  which  should  be  considered  as  "changes 
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of  climate" . Still,  the  geological  record  seems  to  indicate 
that  the  transitions  between  considerably  colder  periods  (ice 
ages  or  glacials)  and  warmer  periods  ("normal"  climates  or 
interglacials)  occurred  over  time  spans  about  ten  times 
shorter  than  the  duration  of  the  relatively  steady  cold  or 
warm  weather  respectively.  This  suggests  what  we  shall  adopt 
here  as  our  operative  definition  of  climate,  viz.,  the  preva- 
lence of  either  warm  weather  (as  we  experience  it  today)  or 
of  cold  weather  (to  mean  a difference  of  the  order  of  ten 
degrees  centigrade  in  yearly  average  below  the  one  recorded 
in  the  present) . 

We  turn  now  to  the  question  of  which  features  of 
atmospheric  phenomena  are  "most  important" . Certainly  tempera- 
ture is  one  of  them,  not  only  because  its  changes  left  deep 
traces  in  the  geological  record  (glaciations  in  temperate 
zones,  pluviations  in  the  tropics  — SMIC,  1971),  but  also 
because  it  affects  all  conditions  of  life  and  because  it  is 
directly  linked  to  the  major  thermodynamic  and  dynamic 
processes  in  the  atmosphere  which  determine  climate  and 
its  changes.  Also,  humidity,  wind  direction  and  intensity, 
cloud  amount,  precipitation,  all  play  a major 
role  in  determining  what  is  perceived  as  weather  and  hence 
should  be  time-averaged  (and,  possibly,  space-averaged) 
into  climate.  Moreover,  it  is  not  only  the  averages  of 
these  quantities,  but  their  day-to-night  and  season-to-season 
contrast  that  enters  our  intuitive  concept  of  a climate. 
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Thus,  at  least  their  variance  should  be  included  in  a more 
complete  mathematical  model  for  climatology. 

In  this  article  we  'hall  treat  a very  simple  model, 
based  on  the  work  of  Sellers  (1969)  and  of  Schneider 
and  Gal-Chen  (1972);  we  hope  that  the  results  will  in 
themselves  be  of  some  significance  for  climate  theory,  as 
well  as  providing  insights  for  devising  and  analyzing  more 
complex  models. 

In  Section  2 the  model  to  be  studied  is  described; 
the  physical  principles  on  which  it  is  based,  as  well  as 
the  empirical  data  it  uses  are  discussed. 

In  Section  3 we  discuss  the  work  of  different  authors 
on  similar  models;  the  similarities  and  differences  between 
their  results  are  pointed  out  and  the  issues  arising  from 
these  results  are  outlined. 

In  Section  4 we  compute  numerically  the  model's 
steady-state  solutions  of  physical  interest,  i.e.,  those 
yielding  positive  absolute  temperatures.  Three  such 
solutions,  corresponding  to  three  distinct  climates  of  our 
planet,  are  obtained:  one  corresponds  to  the  current  climate, 

the  second  to  an  ice  age,  the  third  to  a completely  ice- 
covered  earth.  In  this  section  we  also  explore  the  effect 
of  certain  modifications  in  the  model  on  the  number  of 
steady-state  solutions. 

In  Section  5 the  notion  of  stability  for  the  solutions 
obtained  in  Section  4 is  defined  precisely;  it  is  investigated 
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using  a combination  of  analytical  and  numerical  techniques. 
The  results  are  that  the  present  climate  and  the  ice -covered 
earth  are  stable,  whereas  the  ice  age  of  the  model  is 
unstable . 

In  Section  6 the  effect  of  changes  in  the  solar  radia- 
tion on  the  number  and  stability  of  steady-state  solutions 
is  studied.  The  main  result  is  that  for  a sufficient 
decrease  in  the  solar  radiation  (about  2%) , the  glacial  and 
the  interglacial  solutions  disappear,  leaving  the  ice- 
covered  earth  as  the  only  possible  climate. 


2.  The  Model 


The  model  chosen  for  study  is  a zonally-and-vertically 
averaged  energy-balance  climate  model.  This  means  that 
quantities  in  the  model  are  averaged  over  longitude  and 
height,  leaving  colatitude  <J>  as  the  only  space  variable. 
The  term  energy  balance  means  that  the  model  is  essentially 
based  on  the  energy  equation  of  fluid  dynamics  and  has  sea- 
level  temperature  u as  the  only  dependent  variable.  The 
equation  governing  the  model  is 

(la)  C(4>)ut  = Bi(<|),u)  - R0U,u)  + D(<fr,u,u^,U(|)(J))  ; 

C is  the  heat  capacity  of  the  atmosphere,  land  and  water 
masses;  is  the  heat  absorbed  from  incoming  radiation, 

(lb)  Ri  = Q (4>)  [1  - aU,u)  ] , 

where  Q is  high-frequency  solar  radiation  and  a is  the 
reflectivity  (albedo)  of  the  land  and  sea  surface;  RQ  is 
the  heat  lost  in  outgoing  low-frequency  planetary  re- 
radiation reaching  outer  space, 

4 

(lc)  Rq  = c(<(),u)au  ; 

and  D describes  the  redistribution  of  heat  on  the  surface 
of  the  planet  by  conduction  and  convection, 

<ld>  D - hH  t*tsin  * ' k(‘>'u,lu*  • 

The  coefficients  and  forcing  terms  in  this  model 


represent  yearly  averages  of  the  corresponding  quantities  and 
therefore  do  not  depend  explicitly  on  the  time  t.  Averaging 
the  daily  and  seasonal  variations  seems  justified,  since  the 
time  scales  in  which  we  are  interested  in  our  investigation 
are  of  the  order  of  hundreds  and  thousands  of  years.  The  lack 
of  explicit  dependence  on  time  has  the  advantage  that  the 
model  admits,  as  we  shall  see,  steady-state  solutions,  ufc  = 0, 
which  we  define  to  be  its  climates . The  purpose  of  this  work 
is  to  study  analytically  and  numerically  the  number  of  these 
climates  and  their  stability  under  perturbations. 

The  first  model  of  this  type,  in  finite-difference  form 
and  without  time  dependence,  was  developed  by  Sellers  (1969). 
The  differential  formulation  is  due  to  Faegre  (1972);  time 
dependence  was  introduced  by  Dwyer  and  Petersen  (1973)  and, 
independently,  by  Schneider  and  Gal-Chen  (1973).  Dwyer  and 
Petersen  also  gave  the  outline  of  a systematic  derivation  of 
(1)  from  the  energy  equation  of  fluid  dynamics,  mentioning  the 
main  assumptions  involved.  An  even  simpler  model  has  been 
proposed  by  Budyko  (1969):  in  it  the  diffusion  term  D is 
replaced  by  a nondif ferentiated,  linear  term  in  u,  and  the 
albedo  is  a simple  step  function  of  u only;  this  model  was 
also  discussed  very  thoroughly  by  Leith  (1974),  and  by  Held 
and  Suarez  (1974). 

One  of  the  main  features  of  the  model  (la-d)  is  the  form 
of  the  albedo, 

(2a)  a = {b(<J>)  - c1[um  + (u-c2z  (4>) -um)  _]  )c  , 

where  the  meaning  of  the  subscripts  ( ) _ and  { >c  is 
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given  for  a generic  quantity  h by 


h = min  {h,0} 


0.25, 


h < 0.25, 


h , 0. 25  < h < 0. 85, 

0.85,  0. 85  < h , 


The  subscript  c stands  for  cutoff;  the  cutoff  given  by  (2c) 
embodies  the  observed  minimum  and  maximum  values  of  surface 
albedo. 

Snow  and  ice  have  higher  reflectivity  than  bare  ground 
or  water;  since  in  regions  of  lower  yearly  average  tempera- 
ture the  snow  and  ice  cover  persists  for  a longer  fraction 
of  the  year,  at  lower  temperature  the  yearly  average  albedo 
is  higher;  this  is  expressed  in  the  monotonically  decreasing 
dependence  of  a on  u.  Further,  the  plausible  assumption  is 
made  that,  above  a certain  yearly  average  temperature  um  , 
no  snow  or  ice  will  be  present  at  any  time  of  the  year; 
therefore  a is  independent  of  u for  u— C2 z >_  u^  , as  seen 
from  (2a)  and  (2b).  The  term  c2z(<j>)  gives  the  difference 
between  sea-level  temperature  u and  ground  temperature,  u— C2Z« 
A serious  drawback  of  the  model  is  that  it  does  not 
include  the  effect  on  the  albedo  of  clouds,  atmospheric 
turbidity,  relative  humidity,  and  vegetation.  The  optical 
properties  of  these  factors  and  their  relationship  to  surface 
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temperature  are  less  well  known  and  cannot  be  easily  para- 
meterized in  a model  as  simple  as  the  one  at  hand. 

The  factor  c in  the  outgoing  infrared  radiation 
term  Rq  , 

(2d)  c = 1 - m tanh  (c^uS  , 

expresses  empirically  the  "greenhouse  effect",  i.e.,  the 
screening  by  the  atmosphere,  in  particular  by  the  clouds 
in  it,  of  infrared  radiation  from  the  earth,  thus  preventing 
part  of  it  from  reaching  outer  space.  Notice  that  c decreases 
as  u increases;  this  indicates  that  cloud  formation,  and 
hence  the  opacity  of  the  atmosphere  to  low-frequency  radiation, 
increases  with  increasing  temperature. 

The  function  k(<j>,u)  in  (Id)  has  the  form 

c.  -C-/U 

(2e)  k(4>,u)  = k1(<j>)+k2  (tjj)g(u)  , g(u)  - e = f ' (u)  ? 

u 

k i ( 4* ) u^>  is  sensible  heat  flux  in  the  atmosphere  and  in  the 
oceans,  whereas  k2 (<i>) g(u) u^  is  latent  heat  flux  in  the 
atmosphere.  Here  k^(<(>),  k2(<J>)  are  eddy  dif fusivities , 
which  parameterize  convective  transports;  true  conduction 
is  known  to  be  negligible  in  the  atmosphere-ocean  system 
on  the  planetary  scale.  In  the  original  Sellers  model, 
k(<|>,u)  had  the  form 

k(4>,u)  = kg(4>,u)  - v(u^)  • (u+f  (u) ) , 

where  v is,  for  the  purposes  of  modeling  heat  flux,  mean 
meridional  velocity.  The  theoretical  shortcomings  of  the 
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additional  term  v(u+f(u))  were  pointed  out  by  Robinson 
(1971) ; the  practica1  difficulties  in  giving  a good  para- 
meterization are  discussed  by  Sellers  (1973). 

Numerical  studies  of  Schneider  and  Gal-Chen  (1973) 
indicate  that  results  with  the  original  Sellers  model 
(denoted  by  them  as  (S))  were  very  similar  to  those  with 
the  model  adopted  here  (denoted  in  their  work  as  (SV) ) , 
provided  that  the  numerical  values  of  k and  kg  were 
properly  chosen  (see  the  discussion  on  the  determination 
of  coefficients  further  on).  Furthermore,  the  recent 
work  of  Gal-Chen  and  Schneider  (1975)  shows  that  there  are 
theoretical  grounds  on  which  the  formulation  (SV)  with  zero 
meridional  velocity  is  to  be  preferred.  These  considerations 
will  be  touched  upon  later,  in  a different  connection. 

The  constants  , 1 £ j £ 5,  um  , o , m , as  well  as 

the  empirical  functions  C(<{>),  Q(4>),  b ( 4> ) , z(<j>),  k^(<{>)  and 
k2(<t>)  are  made  to  fit  currently  measured  values  of  temperature 
radiation,  elevation,  albedo  and  heat  flux.  The  functions 
C(<)>),  Q(4>)  and  z(<t>)  are  determined  directly  from  measurements 
The  function  b(<f>)  and  the  constant  c-^  in  (2a)  were 
computed  by  Sellers  (1969)  so  as  to  fit  existing  albedo 
measurements  in  the  northern  and  southern  hemisphere. 

The  constants  m and  c^  were  also  computed  by  Sellers,  so 
as  to  fit  empirical  data  on  RQ ; o is  the  Stef an-Boltzmann 
constant.  The  form  of  the  function  g(u)  and  the  constants 
c 4,  Cj.  appearing  in  (2e)  are  derived  from  certain  physical 


-9- 


i 


considerations  having  to  do  with  the  thermodynamics  of  wet 

air  and  from  corresponding  empirical  data  (see  handbook  of 

Meteorology,  1945).  The  functions  k1(4>),  (<J>)  are  computed 

from  measured  data  on  sensible  and  latent  heat  flux,  k.(d))u 

l $ 

and  k2  ( 4> ) g ( u)  u^  , respectively.  These  computations  are 
based  on  the  measured  temperature  distribution,  denoted 
hereafter  by 

u = u(4>)  , 

which  will  be  called  the  data  climate.  Note  that  we  use 
here  the  term  "climate"  only  for  convenience,  instead  of  the 
lengthier  "temperature  distribution";  u(4>)  is  not  necessarily 
a steady -state  solution  of  the  model;  we  return  to  this  point 
in  Section  4. 

| 

The  measured  data  are  available  at  intervals  of  10° 
latitude  and  are  given  in  Table  1.  Since  the  previously 
quoted  authors  used  finite-difference  formulations  with  a 
fixed  10°  grid  (except  Faegre  (1972),  who  used  a 5°  grid), 
these  data  were  sufficient  for  their  numerical  work.  In  our 
numerical  work,  however,  variable  grid  size  was  employed, 
and  the  10°  data  were  accordingly  fitted  by  smooth  functions. 

In  fact,  in  order  to  have  an  additional  check  on  the  well 
posedness  of  the  model  (i.e.,  the  continuous  dependence  of 
the  solutions  on  the  data,  commonly  referred  to  in  the 
meteorological  literature  as  sensitivity) , two  forms  of  curve 
fitting  were  used;  (i)  by  Bernstein  polynomials,  and  (ii)  by 

i 

cubic  splines.  Results  with  the  two  forms  of  curve  fitting 
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were  indeed  very  similar  (see  Table  2). 

In  this  work  we  only  consider  symmetric  solutions  of  (1) ; 
all  data  are  symmetrized  with  respect  to  the  equator, 

<p  = tt/2.  For  such  data  the  appropriate  boundary  conditions 
are 

(3a)  u^O)  = 0 , (3b)  u^  (n/2)  = 0 ; 

in  the  symmetric  case  these  are  equivalent  to  u, (0)=  u,(n)=  0. 

<P  4> 

We  feel  that  the  slight  asymmetry  between  the  northern 
and  southern  hemispheres  could  hardly  have  had  a major 
influence  on  climatic  change.  Indeed,  the  circulations 
of  the  two  hemispheres  are  practically  separated  from  each 
other  by  the  intertropical  convergence  zone,  which  acts  with 
respect  to  our  model  as  an  insulating  wall.  The  approximation 
involved  in  placing  this  wall  at  the  equator  is  no  worse 
than  other  approximations  in  the  model  (see  also  Held  and 
Suarez,  1974) . A further  reason  for  symmetrization  will 
become  apparent  in  the  next  section. 


3.  Previous  Results 


Budyko  (1969)  and  Sellers  (1969)  used  iterative  numeri- 
cal techniques  for  constructing  solutions  of  their  time- 
independent  models.  They  explored  a range  of  values  of  the 
parameters  appearing  in  the  model,  especially  of  the  solar 
radiation  Q,  and  obtained  one  solution  for  each  set  of 
values.  These  solutions  did  not  depend  smoothly  on  the 
parameter  values;  in  several  instances,  small  changes  in  the 
parameters  led  to  large  changes  in  the  solution.  For 
instance , an  increase  or  decrease  of  a few  percent  in  Q 
resulted  in  temperature  changes  leading  to  extensive  melting, 
or  significant  expansicn  of  the  polar  cap,  respectively. 

Faegre  (1972)  obtained  for  a certain  given  set  of  values 

of  the  parameters  five  distinct  solutions  of  his  variant  of  the 

Sellers  time-independent  model.  Two  of  these  were  highly 

asymmetric,  and  disappeared  when  c(<|>,u)  in  (1)  was  taken 

as  constant;  hence  Faegre  considered  these  solutions  to  be 

spurious  and  unphysical.  It  was  the  desire  to  eliminate  a 

priori  such  solutions  that  suggested  the  choice  of  symmetric 

coefficients.  Faegre ' s formulation  of  the  albedo  is  slightly 

different  from  that  of  Sellers,  mainly  in  that  the  minus 

subscript  in  (2a)  (i.e.,  the  cutoff  of  a(<j>,u)  at  u ) was 

m 

missing . 

The  three  symmetric  solutions  of  Faegre  could  be 
described  as  corresponding  to  the  present  climate,  an 


ice-age  climate  (about  15°  C colder  on  the  average  than  the 
previous  one),  and  a completely  ice-covered  earth  (at  an 
average  temperature  of  about  175  K) . This  last  climate 
was  also  the  one  obtained  by  Sellers  when  decreasing  the 
solar  radiation  by  more  than  4%. 

These  results  raised  the  question  of  the  transitivity 
of  the  earth's  climate,  as  formulated  by  Lorenz  (1968,  1970). 
In  Lorenz's  terminology,  a time-dependent  system  of  equations 
is  transitive  if  its  solutions  have  a unique  equilibrium 
statistic,  that  is,  if  all  solutions,  independently  of 
initial  conditions,  have  the  same  infinite  time  average; 
otherwise  the  system  is  intransitive.  Lorenz  pointed  to  tne 
existence  of  certain  transitive  systems  which  possess  a 
property  called  by  him  almost  intransitivity,  i.e.,  that 
of  having  at  least  some  solutions  whose  averages  over  long, 
but  finite,  time  intervals  are  different  — these  solutions 
then  would  alternate  in  time  between  the  different  averages. 
He  raised  the  possibility  that  the  atmospheric  system  is 
almost  intransitive;  in  other  words,  that  the  known  climate 
changes  in  the  past  were  not  necessarily  caused  by  changes 
in  external  conditions  (like  solar  radiation) , but  rather 
were  an  effect  of  the  normal  evolution  of  the  system. 

Schneider  and  Gal-Chen  (1973)  investigated  the  question 
of  transitivity  for  energy-balance  climate  models  by 
formulating  time-dependent  versions  of  the  Budyko  (B) , 
Sellers  (S  and  SV)  and  Faegre  (F)  models.  They  solved 
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numerically  the  initial-value  problem  governing  these  time 
dependent  models  for  a large  range  of  initial  conditions. 

The  models  were  found  to  be  intransitive,  rather  than  almost 
intransitive:  every  solution  tended  as  t -*■  « to  one  of 

two  (or  more,  in  the  case  of  the  (F)  model)  equilibrium 
solutions;  viz.,  the  equilibrium  statistic  of  the  system 
governing  each  model  was  not  unique,  and  no  spontaneous 
transition  from  one  equilibrium  to  another  was  possible. 

The  two  equilibrium  solutions  obtained  for  all  models 
corresponded  to  the  present  climate  and  to  the  previously 
mentioned  completely  ice-covered  earth. 

The  equilibrium  solutions,  at  least  for  the  (S)  and 
(SV)  models,  proved  stable  under  rather  large  perturbations 
in  both  initial  conditions  and  parameters.  That  is,  solu- 
tions which  differed  in  their  initial  conditions  from  one 
of  the  limiting  "equilibria"  by  as  much  as  + 18  K tended 
as  t +»  to  the  "equilibrium"  near  which  they  started. 

Also  changes  of  + 1.5%  in  the  solar  radiation  led  to  limit- 
ing equilibria  which  differed  by  only  a few  degrees  from 
the  unperturbed  ones.  However,  changes  of  more  than  - 18  K 
i-n  initial  conditions  or  — 1.5%  in  solar  radiation  led  from 
the  "present  climate"  to  the  "ice-covered  earth".  The  latter 
seemed  to  be  the  most  stable  climate  in  all  investigations 
mentioned  above. 

Schneider  and  Gal-Chen  did  not  obtain  a limiting  steady- 
state  solution  which  would  correspond  to  a true  ice  age  as 
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recorded  in  the  planet's  history.  Their  results  with  the 
time-dependent  version  of  the  Faegre  model  (F)  were  rather 
different  from  those  with  the  two  versions  of  the  Sellers 
model  (S  and  SV) , especially  with  regard  to  the  stability 
of  steady  states  under  perturbations. 

Contrary  to  the  results  of  Schneider  and  Gal-Chen, 

Dwyer  and  Petersen  (19  73)  , with  a time "dependent  model 
essentially  identical  to  Schneider  and  Gal-Chen1 s (S)  , 
obtained  solely  one  type  of  limiting  steady  state,  that 
corresponding  to  the  "present  climate".  They  used  only  the 
data  climate  u = u(cj>)  as  initial  conditions,  but  varied 
the  solar  radiation  Q.  The  actual  values  of  the  limiting 
equilibrium  depended  of  course  on  the  values  of  Q used, 
but  slight  changes  in  Q yielded  only  a difference  of  a few 
degrees  between  the  average  temperature  of  the  equilibrium 
and  that  of  the  data  climate;  no  "deep  freeze"  equilibrium 
was  obtained. 

These  results  seemed  to  be  interesting  enough  in  order 
to  warrant  further  study  of  energy-balance  climate  models. 

As  indicated  in  the  previous  section  we  chose  to  investigate 
symmetric  solutions  of  the  (SV)  model  of  Schneider  and 
Gal-Chen,  and  of  some  variations  thereof,  including  the  (F) 
model.  We  hope  that  this  study  will  add  as  much  light  and 
as  little  heat  as  possible  to  the  climate  question  (Jackson, 


We  turn  now  to  the  mathematical  theory  of  equation  (1) . 
Introducing  the  new  space  variable  x = 2<p/v,  we  obtain  the 


initial-and-boundary  value  problem  (IBVP) 


’2.2 


C(x)ut  ' (¥>‘  sin(TTx/2)  sin  ^ •[k1(x)+k2(x)g(u)]ux 


(4) 


- au^tl  - m tanh  (c^u^)] 


+ yQ (x) { l-b(x)+c, [u  + (u-c-z (x) -u  ) ] } , 

^ lm  2.  m-jc 


0 < x < 1 , 0 < t , 


(5a)  u (0,t)  = u (1, t)  = 0 , (5b)  u(x,0)  = u(x)  , 

A A 

where  (4)  is  a nonlinear  parabolic  partial  differential 
equation  (PDE) . Here  g(u)  is  given  by  (2e) , y = 1 
(its  significance  will  appear  later,  in  Section  6)  and 

c^  = 0.009,  C2  = 0.0065  deg  m c^=  1.9  x 10  ^ deg 

3 -2 

(6)  c^  = 6. 105x0. 75xexp( 19. 6) xio  dyn  deg  cm  , c^=  5350  deg, 

a = 1.356x10  ^cal  cm  ^sec  ^deg  m = 0.5,  u = 283.16  deg 

m J 

Mesh  data  at  10°  latitude  for  C(x),  k^(x),  k2(x),  Q(x),  b(x), 
z(x)  are  given  in  Table  1.  The  units  of  the  constants  and 
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mesh  data  are  chosen  such  that  the  common  units  of  all  terms 

...  . -2  -1 
m (4)  are  cal  cm  sec 

The  first  step  of  the  investigation  is  to  find  steady- 
state  solutions  of  (4,  5a)  in  the  range  of  physical  interest, 
and  to  study  their  dependence  on  changes  in  the  model.  We 
consider  therefore  the  steady-state  equation  obtained  from  (4) 
by  setting  ufc  = 0.  After  some  rearrangement  we  get  the 
following  two-point  boundary-value  problem  (BVP)  for  the 
system  of  ordinary  differential  equations  (ODE) : 

(7a) 

(7b) 


(8)  v (0)  = v(l)  = 0 , 

where 

(7c)  k (x,  u)  = k1(x)  + k2(x)g(u)  , 

( 7d)  F (x,  u)  = yQ(x)jl  - b(x)+  c1(um+  (u-c2z  (x)  -uj  _ ] jc 

- au^[l-m  tanh  (c^u^)]  . 

We  use  shooting  (Isaacson  and  Keller,  1966,  Keller,  1968) 
as  the  numerical  procedure  to  solve  (7,  8):  equations  (7a, b) 


u = v 


v = - 


/Ha  2 F(x*u) 
' 2 k(x,u) 


J (cot 


7TX\ 
2 > 


v - 


k^ (x) +k2 (x) g (u) 


k2  (x)  g’  (u)  2 


k (x,  a) 


v 


k (x,u) 


0 < x < 1 , 


v 


are  solved  with  Initial  conditions 

(9a)  v (0)  = 0 , (9b)  u(0)  = uQ  , 

for  different  values  of  Uq;  denote  by  u(x;Ug) , v(x;Uq) 
the  solution  of  the  initial-value  problem  (IVP)  (7,9)  in 
order  to  emphasize  its  dependence  on  the  parameter  Uq. 

An  iterative  scheme  is  then  used  to  obtain  those  values 
of  Ug  which  satisfy 

(10)  v(l;uQ)  = 0 . 

For  these  values  of  Ug  the  solution  u(x;ug) , v(x;Ug) 
of  the  IVP  (7,9)  is  also  a solution  of  the  BVP  (7,8). 

To  obtain  numerical  solutions  of  prescribed  accuracy  to 
the  BVP  (7,8)  one  has  therefore  (Keller,  1968)  to  achieve 
the  desired  accuracy  both  in 
(i)  solving  the  IVP  (7,9),  and  in 

(11)  solving  iteratively  the  nonlinear  (finite)  equation  (10). 
In  solving  numerically  the  IVP  (7,9),  we  encounter  the 

difficulty  that  (7b)  is  singular  at  the  origin,  since 
cot ( ttx/2  ) + « as  x 0.  This  singularity  arises  from  the 
mathematical  form  of  the  diffusion  term  D of  (1)  in  spherical 
coordinates.  It  is  easily  shown,  though,  that  (analytical) 
solutions  of  the  IVP  exist  and  are  unique  at  least  "in  the 
small",  i.e.,  for  | Ug-Ug | < e , 0 < x < 5,  arbitrary  UQ 
and  some  <5,  e >0.  The  numerical  difficulty  of  the  initial 
point  being  singular  can  be  circumvented  by  using  a variable- 


step  method. 


It  also  turns  out  that  in  order  to  achieve  prescribed 
over-all  accuracy  of  the  numerical  solution  of  (7/9)  with 
a minimum  of  computational  effort  it  is  convenient  to  use  a 
variable-step  variable-order  multistep  scheme.  The  ODE 
solver  we  computed  with  was  developed  at  the  Lawrence  Liver- 
more Laboratory  and  is  documented  by  Hindmarsh  (1972,1974). 
Both  the  Gear  and  Adams  methods  included  in  the  package  were 
used  and  gave  results  identical  to  within  the  prescribed 
relative  accuracy,  which  was  10  ; however,  the  Adams  method 

was  faster  and  was  therefore  used  in  most  integrations. 

The  number  of  steps  needed  per  solution  (per  value  of  Uq) 
for  this  accuracy  was  of  the  order  of  100.  It  is  of  interest 
to  point  out  that  near  \_he  singularity  at  x = 0,  and  near 
the  point  where  the  jump  discontinuity  in  the  derivative  of 
F(x,u(x;uq))  occurred,  the  order  of  accuracy  chosen  by  the 
scheme  was  one  and  the  step  size  was  very  small,  so  that  a 
proportionately  larger  amount  of  computation  was  done  in  the 
neighborhood  of  these  two  points. 

Thus  every  solution  of  (7,9)  for  given  u^  puts  at 
our  disposal  a value  of  the  function  v(1;Uq)  , accurate  to 
10-/*.*  To  find  the  zeros  of  v(1;Uq),  we  used  the  me  t hod,  of 
false  position  (Isaacson  and  Keller,  1966).  The  criterion 
for  uQ  to  be  a root  of  (10)  was  |v(1;Uq)  | < 10  . 

* In  the  units  chosen,  typical  values  of  u are  0(10  ) and 
of  v are  0 (10^) . 
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The  range  of  physical  interest  in  which  we  searched  for 


solutions  of  (10)  was  100K  < uQ  < 300K  (see  also  Faegre, 
1972).  Within  this  range,  three  solutions  of  the  BVP  (7,8) 
were  found  which  are  denoted  by  u^(x),  u2 (x) , u^(x). 

They  correspond  to 


u1(0)  = 247.  74  K,  u2(0)  = 223. 97K  , u3  (0)  = 168.94  K. 


The  curve  uxd)  vs.  u(0),  the  zeros  of  which  correspond 
to  the  solutions  u3,u2,u1  , is  given  in  Figure  2.  The 
individual  solutions  u^x),  u2(x),  u3(x)  are  plotted  in 
Figure  3. 


It  is  customary  and  useful  to  characterize  a climate 
of  the  model,  i.e.,  one  of  the  solutions  above,  by  its 
average  temperature,  rather  than  by  the  temperature  at  the 
pole.  Therefore  we  introduce  for  functions  <j>  the 
averaging  operator  6 by 


& (p 


2 /2 
| 4>(x)  sin  (^j)  dx/  j sin  (—)  dx  . 


0 '0 

In  particular,  it  is  known  (Isaacson  and  Keller,  1966)  that 
for  functions  <p  symmetric  about  the  equator,  <J>(2-x)  = <}>(x)  , 
which  also  satisfy  =0,  and  hence  can  be  extended  so 

as  to  have  period  2,  the  trapezoidal-quadrature  approximation 
of  & is  very  accurate.  Denote  this  approximation  by  , 
where  A = 1/9  (corresponding  to  a 10°  mesh)  ; then  we  have 


S.u.  = 287.76  K,  &.u-  = 267.44  K,  &.u,  = 175.43  K . 


Clearly  u^  corresponds  to  an  interglacial,  u2  to  a glacial, 


and  to  the  ice-covered  earth. 


Note  that  for  the  data  climate  u = u(x) 


u (0)  = 247.36  K,  &Au  = 287.20  K , 


which  is  indeed  very  close  to  the  interglacial  solution  of 


the  model,  (see  also  Figure  3 > . Though  present  observa- 


tions were  used  to  obtain  constants  and  empirical  functions 


in  (4) , no  explicit  term  was  added  to  ensure  that  (4)  hold 


with  ufc  = 0?  also  the  diffusion  term,  D,  is  of  the  same 


order  of  magnitude  as  the  radiation  terms,  and  , so 


that  this  is  not  the  result  of  a simple  radiation  balance. 


The  same  is  true  of  the  work  of  all  the  previous  authors 


discussed  in  Section  3 ? these  authors,  however,  assumed 


at  least  implicitly  that  the  data  climate  should  be  a 


steady-state  solution  of  the  model,  or  approximately  so. 


Since  this  result  is  not  actually  built  into  the  model,  as 


far  as  we  can  see,  we  think  it  is  rather  remarkable  that 


this  class  of  models  yields  a steady-state  solution,  u^(x), 


very  close  to  the  data  climate,  u(x).  Henceforth  we  shall 


refer  to  u^  also  as  the  "present  climate"  of  the  model. 


In  fact,  Schneider  and  Gal-Chen  (1973)  did  use  an  extra 
"fudge  coefficient",  0.97  £ cQ(x)  i 1.03,  in  RQ  in  order 
to  achieve  better  agreement  between  the  interglacial,  or 
"present"  climate  of  their  models  and  the  data  climate; 
they  obtained  £u^  = 287.06  K vs.  £u  = 287.30  K for  the 
( S ) mode 1 . 
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The  mesh  data,  from  which  constants  and  empirical 
functions  for  the  model  were  computed,  are  known  to  be  in 
error  by  possibly  as  much  as  100%  (Schneider  and  Gal-Chen, 
1973);  also  some  of  the  parameterizations  used  in  formulati  ng 
the  model  are  questionable.  Therefore  we  made  a number  of 
computations  in  order  to  obtain  information  on  the  dependence 
of  the  physically  significant  solutions,  u^u^u^  on  varia- 
tions in  the  model.  The  salient  features  of  these  computa- 
tions, summarized  in  Table  2,  are: 

(1)  The  dependence  of  all  these  solutions  on  the 
coefficients  seems  to  be  very  smooth. 

(2)  The  bounds  on  the  albedo,  0.25  £ a(x,u)  £0.85,  are 
essential  for  the  existence  of  three  steady-state  solutions 
in  the  physical  range,  100  K £ u(0)  £ 300  K:  u3  disappears 
when  the  bound  a £ 0.85  is  not  enforced,  and  u^  disappears 
when  the  bound  0.25  £ a is  not  enforced. 

(3)  The  terrestrial  radiation  term  Rq  , given  by  (lc), 

and  its  nonlinearity  are  essential  for  the  existence  of 
physically  significant  solutions:  if  the  term  is  set  equal 

to  zero  all  solutions,  , disappear;  when  it  is 

replaced  by  its  average, 


1 

1 


i 


i 


} 

A 


~ 4 

Rq  = £^c(x,u(x))au  (x)  = const.  , 

u2  only  obtains;  when  Rq  is  replaced  by  a linear  approximation, 
Rq  = au^{u+4(u-u) }£^c(x,u(x) ) , u = £^u(x)  , 
then  u^  only  obtains. 
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(4)  Our  results  with  the  model  using  the  Faegre  formula- 
tion of  the  albedo  are  very  similar  to  those  using  the  Sellers 
formulation;  it  is  hard  to  explain  the  disagreement  in  this 
respect  between  our  results  and  those  of  Schneider  and 
Gal-Chen  (1973),  who,  as  mentioned  in  Section  3,  obtained  very 
different  results  when  using  the  two  albedo  formulations. 

(5)  The  values  of  the  empirical  functions  in  the  model 
when  fitting  mesh  data  by  (i)  Bernstein  polynomial  approxima- 
tion and  <ii)  cubic  spline  interpolation  are  pointwise  rather 


I 


different  for  some  of  the  functions  (see  Figure  1) . The 
solutions  of  (7,8)  when  using  these  different  fits  are,  however, 
very  close.  This  also  supports  the  assertion  in  paragraph  (1) 
above,  and  shows  that  the  uncertainty  in  the  mesh  data  does 
not  affect  in  an  important  way  the  conclusions  of  the  investi- 
gation. 

We  explored  another  variant  of  the  model,  in  which  the 
eddy  diffusivity  k(x,u)  given  in  (7c),  corresponding  to  the 
(SV)  model  of  Schneider  and  Gal-Chen,  was  replaced  by  k(x,u), 
where  u is  the  data  climate;  hence  also  g' (u)v  in  (7b)  is 
replaced  by  g'(u)u'(x).  This  variant,  in  accordance  with  the 
terminology  of  Gal-Chen  and  Schneider  (1975),  would  correspond 
to  an  (SVC)  model.  We  denote  this  modification  of  (7.)  by  (7'); 
the  solutions  of  (7'  ,8)  are 

u.  (0)  = 247.55  K,  u-(0)  = 227.76  K,  ^(0)  = 169.44  K, 

&Aux  = 287.  70  K,  £au2  = 268.60  K,  £^3  = 175.44  K. 
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These  are  indeed  very  close  to  the  solutions  of  (7,8), 
especially  for  and  (see  also  Figure  3) . 

Gal-Chen  and  Schneider  (1975)  investigated  the  effect  of 
variations  in  the  solar  radiation  yQ  on  the  equator-to-pole 
temperature  gradient;  they  argue  that  this  dependence  should  be 
monotone,  in  fact  monotone  increasing.  In  the  light  of  this 
argument  and  of  their  results,  model  (SV)  is  superior  to  (S) , 
as  already  mentioned  in  Section  z;  moreover,  model  (SVC)  is 
superior  to  (SV) . Since  (SVC)  is  also  more  convenient  to  use 
when  investigating  the  stability  of  the  steady-state  solutions, 
we  shall  work  with  it  in  the  sequel;  the  corresponding  form 
of  (4)  we  denote  by  (4')/  the  corresponding  form  of  (1)  by  (1') 
and  the  corresponding  form  of  (7)  by  (7'). 

For  this  model,  additional  computations  were  made  outside 
the  range  of  physical  interest.  One  more  steady-state  solution 
was  found;  we  denote  it  by  u^(x),  and  it  satisfies 

u4(0)  = - 185.99  K,  £au4  = - 175.40  K . 

Most  probably  there  are  no  other  solutions  of  the  BVP  (7'  ,8) 
at  all.  Indeed,  the  solutions  of  the  IVP  (7', 9)  become 
unbounded  as  Uq  moves  towards  the  ends  of  the  range  explored, 
viz.,  -1330  K < uQ  £ 300  K;  i.e.,  u(x;uQ)  -*•  +»  as  uQ  approaches 
the  ends  of  the  interval  above  (see  Figure  2). 

We  would  like  also  to  point  to  the  fact  that  some  of 
the  modifications  of  the  model  which  have  three  physical 
steady  states  yielded  numerical  values  of  the  latter  consider- 


i ■■  — - — • --  '—■*-*•“  r-  **■“ — 

\ 

\ 

ably  higher  than  (see  Table  2);  however,  it  was  the 

qualitative  feature  of  the  number  of  solutions  and  their 
approximate  position  with  respect  to  each  other  that  was  of 
interest  in  our  investigation:  the  averaae  temperature  of 

a given  solution  of  any  fixed  model  could  be  changed  to 
practically  coincide  with  that  of  the  solution  Uj  of  (4)  to 
which  it  corresponds  by  adjusting  the  values  of  the  coeffi- 


i 


cients . 


5.  Stability  of  the  Steady-State  Solutions 


In  the  previous  section  we  have  shown  that  equation  (4') 
with  the  boundary  conditions  (5a)  has  three  steady-state 
solutions  of  physical  interest,  u = u^  (x)  , 1 ^ j _<  3.  In 
this  section  we  shall  study  the  stability  of  these  steady 
states . 

Let  us  concentrate  on  any  one  of  the  steady  states  above, 

u = Uj (x) , where  j is  fixed.  Stability  of  u^  means  that  small 

perturbations  of  u^  die  out  with  time.  More  precisely,  u. 

0 

is  stable  if,  when  taking  any  nearby  state  U,  i.e. , one 

which  differs  little  from  u . , 

3 

O O 

( 11)  U = Uj  (x)  + £ V ( x)  , 

O 

say,  where  v is  arbitrary  and  c > 0 is  not  too  large,  then 
the  solution  U(x,t;e)  of  (4' -5)  with  initial  condition 

O 

U(x,0;e)  = U(x;e)  tends  to  u^  itself  as  t -*■  00  . 

Let  (4')  be  written  symbolically  as 

(12)  ufc  = N (u)  , 

where  we  divided  through  equation  (4')  by  C (x)  , 0 < Cm  _<  C(x) 

_<  CM  < ® , and  N(u)  is  the  corresponding  right-hand  side.  Take  a 
perturbed  u,  u = U(x,t,e),  which  is  assumed  to  satisfy 
the  boundary  conditions  (5a)  and  the  initial  condition  (11). 

If  such  a solution  of  (12)  exists  for  all  e sufficiently  small, 
0 < e £ £q  , then  the  equations 


(12*) 


Ufc (x,t;e) -Ut (x, t ;0)  N (U  (x,  t ; e ) ) -N  (U  (x,  t; 0 ) ) 


0 < e <_  eQ  , 


also  hold,  with  U(x,t;0)  = u^ (x) . Letting  e -*•  0 we 
obtain  the  linearized  equation 


(13a)  v,  = - L .v  , 

^ D 

j 

where  we  define  v(x,t)  = U(x, t;e) | e_q  and 


(13b)  Lj  = ‘ N{u(x,t;e)) |e=0  = - f^Nlull  , l<j<3. 

At  this  point,  for  the  sake  of  simplicity,  we  shall 
drop  the  subscript  j,  and  u will  thus  stand  for  some  u^  , 
L for  the  corresponding  L ^ • 

Equation  (13a)  is  linear  in  v and  it  has  a unique  solu- 
tion v satisfying  the  boundary  conditions 

(13c)  v (0,t)  = v (l,t)  = 0 , 

X X 

and  arbitrary  initial  conditions 
(13d)  v(x,0)  = v(x)  . 

The  solution  may  be  found  by  the  usual  method  of  separation 
of  variables,  or  expansion  in  eigenfunctions  (normal  modes). 
Consider 

v(x,t)  = e ^ w(x)  ; 
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this  v will  be  a solution  of  (13)  with  v(x,0)  = w(x)  iff 
(if  and  only  if)  A is  an  eigenvalue  and  w is  an  eigenfunction 
of  L,  i.e.,  iff  the  pair  (A,w)  satisfies  the  homogeneous 
problem 

(14a)  Lw  = Aw 

(14b)  w (0)  = wx(l)  = 0 . 

We  shall  show  that  the  operator  L has  sufficiently  many 
eigenfunctions,  and  that  therefore  any  solution  v of  (13) 
can  be  written  as  a series 

r -a  ^)  t (if) 

(15)  v(x,t)  = l a^e  z wlKJ  (x) 

k=0  K 

(k)  (k) 

where  (A  ,w  ')  are  all  the  solutions  of  (14).  It  is 

clear  from  the  derivation  (12')  of  (13)  and  from  (15)  that 

for  the  steady  state  u to  be  stable  it  is  necessary  that 

(k) 

every  eigenvalue  A of  L have  positive  real  part: 

(16)  Re  A(k)  > 0 . 

This  condition  defines  the  so-called  linear  stability  of  u. 

It  has  been  shown  rigorously  in  a few  cases  and  it  is  believed 
in  many  others  that  (16)  is  also  sufficient  for  the  stability 
of  u as  we  defined  it  at  the  beginning  of  this  section,  i.e., 
that  linear  stability  implies  nonlinear  stability.  In  the 
next  section  we  shall  give  an  argument  which  will  make  it  at 
least  plausible  that  this  is  the  case  also  for  the  problem 
at  hand. 
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We  turn  now  to  the  analysis  of  the  eigenvalues  of  the 
linear  second-order  ordinary  differential  operator  L,  which 
we  can  write  as 


(17a) 

Lw  = - 

r(x)  (P(x)wx)x  + <3(x)w  ‘ 

Here  p 

,q,r  are  determined  by  (13b)  and  by  (4'): 

(17b) 

r(x)  - C(x)  sin  ( 2 ) , 

(17c) 

p(x) 

= (ijj-)  2 sin  (— j^)*k(x,u)  , 

(17d) 

q (x)  = 

(Q(x)(c1)c,  - d (x,  u)  }/C  (x)  , 

wi  th 

f c.^  if  0.25  £ a(x,u)  £ 0.85 

( 17e) 

(cl>c'  * ' 

and  u(x) -c2z  (x) -u^  £ 0 , 

0 otherwise, 

and 

( 17f ) 

3 . . 4 

d(x,  u)  = c(x,u)  au 

= ou^|4[l-m  tanh(c2u^)]-  6mc2u^  [l-tanh^  (c^u^)  ] j 

Clearly  L thus  defined  is  formally  self-adjoint  (Courant 
and  Hilbert,  1953,  Birkhoff  and  Rota,  1969)  under  the 
prescribed  boundary  conditions  with  respect  to  the  inner 
product 

(18)  (w1,w2)  = | r(x)w1(x)w2  (x)  dx  . 

0 


Moreover,  r,p  are  nonnegative  and  twice  continuously  differ- 


entiable  on  the  interval  I = [0,1],  and  q is  piecewise  contin- 
uous on  I.  However,  because  of  the  singularity  at  x = 0 due  to 
the  fact  that  r(0)  = p(0)  = 0, the  usual  Sturm-Liouville  theory 
of  self-adjoint  operators  (Courant  and  Hilbert,  1953,  Birkhoff 
and  Rota,  1969)  does  not  apply  to  L;  this  difficulty  though 
can  be  overcome  and  the  theory  can  be  extended.  The  crucial 
element  in  this  extension  is  the  observation  that,  because 
of  the  boundedness  of  q,  L is  bounded  from  below  in  the  sense 
that,  for  any  w satisfying  (14b)  for  which  (w,w)  < °°,  the 
inequality 

(19a)  <w,w>  >_  K(w,w)  , 

holds  with  some  fixed  constant  K,  I<  > min  q(x)  >-<*>,  independent 

0£x<l 

of  w;  here  <*r,w>  is  the  Dirichlet  integral  corresponding  to  L, 

1 

f 2 2 

(19b)  <w,w>  = (Lw,w)  = (Pwx  + rc3w  ) dx  • 

0 

This  result,  together  with  an  analysis  of  the  nature  of 
the  singularity  at  x = 0,  are  sufficient  to  show  that  indeed 
L has  a complete  system  of  eigenfunctions,  orthonormal  with 
respect  to  the  inner  product  (18)  (Birkhoff  and  Rota,  1969)  , 
and  that  the  eigenvalues  of  L are  real  and  can  be  arranged 
in  an  ascending  sequence 

- < A(1>  < A<21  < A(3>  < ... 

with  A(k)  -*■  00  as  k -*■  00 , and  K = A^  in  (19).  Hence 
any  solution  v of  (13)  can  be  written  as  a series  (15) . 
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Therefore  the  stability  question  for  the  steady  state  u 
reduces  to  determining  whether 

(16')  X(1)  > 0 , 

in  which  case  u is  stable,  or  whether  the  opposite  holds 
in  which  case  u is  unstable . 


5a.  Stability  Criterion 


In  this  subsection  we  shall  show  that  it  is  possible  to 
determine  the  sign  of  A^^  without  actually  computing  A^. 

We  start  with  the  variational  characterization  of  the 
lowest  eigenvalue  (Courant  and  Hilbert,  1953),  A(1), 

, (1)  . <v,v> 

(20)  A = min  r-  = min  <v,v>  , 

<v,v><°°  v,v  ( v , v)  =1 

( v , v ) 

i.e.,  A^  is  the  minimum  of  the  Rayleigh  quotient  R(v) 
corresponding  to  L, 

(21)  R(v)  = <v,v>/(v,v)  , 

and  the  minimum  is  assumed  for  v = w^.  Indeed,  (14a)  is 
the  Euler  equation  for  (20),  and  vx(l)  = 0 is  the  "natural" 
boundary  condition  for  (20)  in  the  sense  of  the  calculus  of 
variations  (Courant  and  Hilbert,  1953) ; also  v (0)  = 0 is, 
according  to  the  theory  of  singular  operators  with  the 
properties  of  L,  the  only  boundary  condition  at  x = 0 which 
ensures  that,  for  solutions  v of  (14),  <v,v>  as  well  as  (v,v) 
is  finite.  In  particular,  we  conclude  that  the  minimizing 
function,  v = w^,  is  at  least  as  smooth  as  the  coefficients 
of  L,  viz. , it  must  have  a piecewise  continuous  second 
derivative . 

With  these  preliminaries  we  are  able  to  prove  the 
following  known 
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Lemma. 


The  first  eigenfunction  of  L,  w 


, is  strictly 


(1) 


positive , 

w ^ (x)  > 0 , 0 x £ 1 . 

Proof ; From  the  definition  (21)  of  R(v)  it  is  clear  that 

R(  |w(1)  | ) = R(w(1)  ) , 

where  |y|  denotes  the  absolute  value  of  y.  If  w'1'  had  a 
zero  at  some  interior  point  Xq  , 0 < Xq  < 1,  and  w^^  (Xq)^O, 

then  the  first  derivative  of  |w^|  would  have  a jump  at 
x = Xg  , which  contradicts  the  smoothness  of  |w^|  as  a 
solution  of  the  variational  problem.  If,  on  the  other  hand, 
w“>  (xQ)  =0,  0 < x0  < 1,  or  w(1)  (0)  = 0,  or  w(1)  (1)  = 0, 
then,  by  che  uniqueness  theory  of  linear  ODE  (Birkhoff  and 
Rota,  1969),  it  would  follow  that  w^  5 0;  this  contradicts 
the  fact  that  w^  is  a nontrivial  solution  of  (14), 
completing  our  proof. 

Now  we  are  ready  to  state  our  stability  criterion,  which 
is  part  of  the  conclusion  of  the 

Theorem.  Let  L,  \ = X^,  w = w^  be  as  above.  Suppose 
also  that  there  exists  a function  v,  as  smooth  as  w, 
nonnegative,  v >_  0 on  I , and  satisfying 

(22)  Lv  > 0 , v (0)  = 0 , v (0 ) = 1 . 

X 

Then  v (1)  >_  0 implies  X _>  0 . Moreover, 

X 

(i)  X > 0 if  either  Lv  % 0 or  v (1)  > 0,  and 

X 

(ii)  if  Lv  = 0 holds,  then  X>0,  A=0,  or  A<0,  according  to 

whether  v (1)  >0,  v (1)  = 0 , or  v (1)  < 0. 
Xx  x 
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Proof:  The  required  results  are  easily  read  off  from  the 

following  sequence  of  eq  .ities  obtained  from  the  definitions 
and  by  integration  by  parts: 


X j rvw  = 
0 


A(w,v)  = (Lw,v)  = ~(PWx^xV  + r<3wv 


1 1 


-pw  V 
X 


= pw  V 


+ J pw^v^  + rqwv 


0 0 
1 


r 


+ 

0 0 


"(pv  ) w + rqvw 

A A 


= (pwvx) (1)  + (Lv,w)  , 


where  we  used  w (0)  = w (1)  = 0 and  v (0)  = 0 . 

XX  X 

Note.  This  is  essentially  a comparison  theorem,  of  the 
type  familiar  from  S turm-Liouvi lie  theory  (Birkhoff  and  Rota 
1969)  . 


The  result  under  (ii)  above  is  the  stability  criterion 
which  we  shall  use.  For  completeness  we  give  here  also  the 
following  slight  generalization  as  a 


Corollary.  Let  L,  X,  w be  as  in  the  Theorem.  Suppose  v 
satisfies  the  hypotheses  of  the  theorem,  except  for  that  of 
being  nonnegative  on  I.  Instead,  let  x^  , 0 < x^  <_  1,  be 
the  first  zero  of  v, 

v(x)  > 0 on  0 <_  x < Xq. 

Then  X < 0 . 
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Proof:  The  result  follows  from  the  same  integration  by 

parts  carried  out  in  the  proof  of  the  theorem,  except  that 
now  the  upper  limit  of  integration  has  to  be  x^  rather  than 
1,  and  we  use  v(xn)  = 0 rather  than  w (1)  = 0 in  passing 
from  the  second  to  the  third  line.  Moreover,  since  v(x)  > 0 
for  x < Xq  , and  v is  continuously  differentiable,  v^(Xq)£  0. 
Furthermore,  it  cannot  be  that  both  Lv  = 0 and  v^(Xq)  = 0, 
since  then  we  would  have,  by  uniqueness,  v = 0,  which 
contradicts  v(0)  = 1.  This  completes  the  proof. 


5b . Stability  Results 


In  this  subsection  we  shall  apply  the  stability  criterion 
obtained  in  Subsection  5a  to  the  steady -state  solutions  u.  , 

1 < j < 3. 

For  this  purpose  we  construct  functions  Vj  by  solving 

(22)  with  the  equality  sign,  and  with  L « L_.  given  by  (13b) 

and  (17).  The  results,  obtained  with  a relative  numerical 
- 7 

accuracy  of  10  , are  that 


(a) 

Vj  , j = 1,2,3,  is  positive, 

^ > V . > 0 , Vj 

1 v 
o 
Ln 

■*  a 

(b) 

vj (1)  > 0 for  j = 1,3, 

whereas  v.(l)  < 

dx  j 

0 for  j=2 

The 

actual  computed  values  are 

^ v^l)  = 2.39970,  ^ v2(l)  = -3.24312,  (1)  = 4. 31025 . 

These  values  of  the  derivatives  at  x = 1,  together  with  the 
values  of  Vj , are  sufficiently  bounded  away  from  zero  in  order 
to  conclude  that  the  conditions  of  the  theorem  are  satisfied 
beyond  the  doubt  of  numerical  uncertainty. 

Hence  the  solutions  representing  the  present  climate  and 
the  ice-covered  earth  are  stable,  whereas  the  so-called  ice 
age  of  the  model  is  unstable.  This  result  agrees  with  and 
throws  additional  light  on  the  results  of  previous  authors,  in 
particular  the  time -dependent  integrations  of  Schneider  and 
Gal-Chen  (1973). 
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Actually  the  stability  computations  were  carried  out 
also  for  the  solutions  of 


N 1 (Uj ) =0  , 


1 1 j 1 3, 


with 


(23)  N'  (u)  = £)2  — : — . — — i-  sin  (2-Z)  -u 

tt  sm  (ttx/2  ) 9x  ' 2 1 x 


Q(xl{l  - Bq+  clU}0  - 


c,au 

6 


where  the  subscript 


is  defined  in  (2c)  and 


K0  = SAk<x»u(x))  = 2.2*10  , BQ=  £Ab(x)  = 2.85881, 


c6  = fiAc(x,u(x))  = 0.61  ; 


tile  corresponding  linearization  of  N'  is 


(13')  L'  = - N (u.)  = (V  ,°  . 1-  sin  f— ) 

3 u j tt  sin  ( ttx/2 ) 9x  2J  3x 


where 


+ Q(x)  (c. ) ,,  - 4c, au.  , 
1 c 6 j 


<cl'c"  = 


c1  if  0.25  _<  1-BQ+C1u^  _<  0.85, 


0 otherwise 


This  is,  according  to  the  experiments  carried  out  in  Section  4 
(see  Table  2),  the  simplest  form  of  (4)  which  still  yields 
approximately  the  same  steady -state  solutions  in  the  physical 
range  of  interest.  The  results  are  the  same,  i.e.,  u1  and  u3 
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are  stable,  and  is  unstable.  It  seems  therefore  quite 
plausible  to  conclude  also  that  for  all  models,  lying  in  some 
sense  between  (4)  and  (4")  , the  steady  states  corresponding 
roughly  to  ui'u2,u3  have  the  same  stability  properties  as  the 
latter.  Combining  this  remark  with  the  one  made  at  the  end  of 
Section  4,  it  seems  that  we  have  a result  about  the  stability 
of  the  steady  states  for  a certain  type  of  energy-balance  model, 
rather  than  just  for  one  specific  moiel  of  this  type.  It  seems 
desirable  to  define  precisely  the  type  of  model  having  these 


properties,  and  we  intend  to  try  to  do  so  in  further  work. 

Having  thus  determined  the  stability  properties  of  the 
steady  states  of  (4'),  we  want  to  obtain  some  additional  infor- 
mation by  actually  computing  the  lowest  eigenvalues  xl1*,  and 
• ( 1 ) 

corresponding  eigenfunctions  w_j  , of  L ^ , 1 <_  j <3.  The  eigen- 
values are 


| 

i 


x{1}=  3.95390xio"9,  X^1}=  -5. 87662><lo"9  , = 5 . 13587x  lo'9  , 

and  the  eigenfunctions  wj1*,  1 < j < ^ are  plotted  in  Figure  4. 

The  method  for  computing  (xj1*,wj1*)  was  again  shooting, 
this  time  with  respect  to  the  parameter  A.  That  is,  equation  (14a) 
with  L = Lj  , was  solved  for  w(x;A)  with  initial  conditions 
w ( 0 ; A ) = 1 , w ( 0 ; A ) = 0 , 

and  with  different  values  of  the  'shooting  parameter"  A.  The  zeros 
of  the  function  vx(l;A)  were  then  computed  by  regula  falsi  to 
within  an  accuracy  of  10  . The  over-all  relative  error  in 

computing  solutions  of  the  initial-value  problem  was  10_^, 
as  before. 


3 


I 

i 

: 
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Notice  from  (15)  that 

x = 1/1  X j1)  | 

is  a characteristic  decay  or  relaxation  time  for  the  solution 
U(x,t;e)  of  (12)  to  u = Uj(x).  This  time  is  of  the  order 
of  10  years  for  all  j, 

x1  = 8.0  years  , x2  = 5.4  years,,  x^  = 6.2  years. 

In  this  context  it  is  remarkable  that  Schneider  and  Gal-Chen 
(1973)  state  that,  in  one  of  their  integrations,  the  solution 

O 

of  (4)  with  initial  data  u(x)  = u(x)  and  y = 0.984,  after 
dropping  rapidly  by  about  12  K in  average  temperature,  was 
nearly  constant  for  about  50  years  of  simulated  time,  and 
then  finally  dropped  to  a steady  state  close  to  our  u^(x). 

From  our  results  it  becomes  clear  that  the  solution  mentioned 
above  hovered  for  a time  of  the  order  of  x2  around  u2  , but 

could  not  persist  there  indefinitely  because  of  the  instability 

* 

of  u2  , and  finally  attained  u^  • which  was  stable. 

It  also  follows  from  (13)  and  (14)  that  multiplying  C (x) 
by  a constant  < > 0 will  result  in  x^,  1 < j £ 3,  being 

multiplied  by  <•  A similar  statement  holds  also  for  the 
nonlinear  equation  (4) , since  such  a constant  k just 
corresponds  to  a different  scaling  of  the  time  t (see  also 
Schneider  and  Gal-Chen,  1973).  Experiments  in  determining 

We  shall  see  in  the  next  section  that  x2  increases  with 
decreasing  y. 
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characteristic  response  times  for  solutions  of  (4')  were 
done  by  Dwyer  and  Petersen  (1973),  who  used  two  heat  capaci- 
ties C(x),  both  larger  than  that  used  by  Schneider  and 
Gal-Chen  (1973)  and  by  us.  It  seems,  however,  that  upon 
decreasing  y in  (4)  to  y = 0.98,  as  they  always  took 

O 

u(x)  = u(x)  in  (5b) , the  average  temperature  of  the  solution 
u(x,t)  dropped  rapidly  at  first,  and  only  slowly  thereafter, 
as  indicated  by  Figure  2 in  their  article.  Apparently  this 
slow  decrease,  which  shows  that  their  solution  u(x,t)  of  (4') 
was  approaching  a steady  state  close  to  our  U2(x),  was 
interpreted  by  them  as  proving  the  nonexistence  of  a steady 
state  u^(x) ' which  had  been  obtained  in  the  work  of  Budyko, 
Sellers  and  Faegre  when  decreasing  y by  a similar  amount. 

The  influence  of  changes  in  y on  the  steady-state 
solutions  Uj  , 1 <_  j £ 3,  of  (4')  will  be  investigated  in 
the  following  section.  At  this  point  we  only  want  to  mention 


that,  whereas  a multiplicative  factor  k in  C(x)  affects 

(1) 


the  magnitude  of  A. 


and  hence  of  , 


our  theorem  shows 


that  it  does  not  affect  the  actual  stability  of  u^  , i.e., 


the  sign  of  A • 


(1) 


6 . Perturbed  Steady-State  Solutions  and  Bifurcation 


It  is  clear  that  steady-state  solutions  of  (4')  could 
not  spontaneously  evolve  into  each  other.  More  precisely, 
the  solution  of  (4'-5a)  with  initial  condition  u(x,0)  = Uj(x), 
j = 1,2,3,  is  u(x,t)  = Uj (x) . Moreover,  it  stands  to 

O 

reason  that,  for  any  physical  initial  condition,  u(x)  100  K, 
we  would  have 

lim  u(x,t)  = u.(x)  , 

t-KX>  J 

with  j = 1 or  j = 3,  since  ^2 (x)  is  (linearly)  unstable, 
whereas  u1(x),  u^tx)  are  (linearly)  stable  (and  u^ (x)  <-170K 
< 0)  . 

Thus,  to  explain  physical  ice  ages,  one  has  to  consider 
perturbations  in  the  parameters  appearing  in  equation  (4'). 
Such  perturbations  would  presumably  be  caueed  by  physical 
mechanisms  not  included  in  the  model.  The  most  likely 
candidate  for  such  a role  is  vi,  which  up  to  now  was  taken 
to  be  unity.  Indeed,  many  ice-age  theories  rely  heavily  on 
a change,  however  small,  in  the  amount  of  solar  radiation 
reaching  the  lower  layers  of  the  atmosphere  (SMIC,  1971, 
Beckinsale,  1965) . 

Some  attribute  the  assumed  decreases  in  solar  radiation 
to  changes  in  the  parameters  of  the  motions  of  our  planet 
(Milankovitch,  1930) , others  to  airborne  volcanic  dust 
due  to  an  increase  in  volcanic  activity  (Fuchs  & Patterson, 
1947),  and  so  on.  There  has  also  been  concern  about  a 
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possible  climatic  catastrophe  being  imminent  because  of  the 
increase  in  the  quantity  of  industrial  pollutants  in  the 
atmosphere  (Rasool  and  Schneider,  1971). 

To  investigate  the  effect  of  such  changes  in  the  model 
at  hand,  the  curve  u (1)  vs.  u(0)  was  recomputed  for 
different  values  of  y,  in  particular  with  a view  to  obtain- 
ing  u1(x;y),  (x; y) . One  important  result  is  that  these 

two  solutions  coallesce  for 

yc  = 0.98152 

and  disappear  entirely  for  y < y . The  bifurcating  solution 

u = uc(x)  = = u2^x;tic^  was  comPute(^  with  over-all 

relative  accuracy  10  and  |uc(l)|  < 5*10  . The 

(u(0) ,u  (1) ) -curve  corresponding  to  y = y is  given  in  Figure  5; 
**  c 

notice  that  it  is  very  flat  near  u(0)  = u (0),  which  makes 

c 

it  difficult  to  compute  uc(x)  accurately. 

Further  computations  were  carried  out  for  y = 0.982, 
0.985,  1.01,  1.02,  1.03,  1.04  and  1.05.  The  results  of  these 
computations  are  summarized  in  Table  3 and  plotted  in  Figure  6. 
It  is  quite  interesting  that,  whereas  for  y _>  1.0  the 
dependence  of  u^(0),  fi^u^ (x) , j = 1,2,  on  y is  almost  linear, 
this  dependence  is  definitely  quadratic  near  y = yc  ; the 
latter  is  in  good  agreement  with  the  general  theory  of 
bifurcation  for  nonlinear  parabolic  problems  (Hoppensteadt 
and  Gordon,  1975).  According  to  the  theory,  there  exists 
in  principle  the  possibility  that,  instead  of  disappearing 
at  y = yc  , the  two  solution  branches  u1(x;y),  U2(x;y)  could 
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merge  into  one  periodic  solution  u^2(x,t;y)  for  y < yc* 
This  possibility  was  not  borne  out,  however,  by  the  time- 
dependent  computations  of  Dwyer  and  Petersen  (1973),  and  of 
Gal-Chen  and  of  Schneider  (1973,  1975);  we  did  not  investi- 
gate it  further. 

Concerning  the  problem  of  the  pole- to- equator  tempera- 
ture gradient, 


AUj  = Uj (1)  - Uj  (0)  , 


1 1 j £ 3, 


as  discussed  by  Stone  (1973)  and  Gal-Chen  and  Schneider  (1975), 
the  curve  Au.  = Au^(y),  j = 1,2,  is  particularly  interesting. 

We  notice  that 

(a)  the  values  of  Au1  lie  below  those  for  Au2  ; 

(b)  the  values  of  Au1  are  monotonically  decreasing  with  y; 

(c)  the  values  of  Au2  have  a maximum  for  y somewhere  between 

y = 0.985  and  y = 0.99; 

(d)  the  dependence  of  both  Au-j^  and  Au2  on  y,  but  especially  that 
of  Au-^  , is  very  steep  near  y = yc* 

Being  aware  of  the  limitations  of  the  model,  as  pointed  out 
in  Section  2,  we  do  not  want  to  make  extensive  comments 
concerning  these  results , but  only  note  them  for  comparison 
with  the  results  of  other  models  and  for  further  study. 

We  also  studied  the  stability  of  the  solutions  u = uc(x) 
and  u = u,(x;y  ).  Repeating  the  construction  indicated  at 
the  beginning  of  Subsection  5b  for  v = vc  ar|d  v = v3  ' where 
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L . = - -r—  N ( u ; y ) I , 

i 3u  1 u=u . 

3 

we  obtained  the  following  results: 


j = c,3. 


(a)  v.  , j = c,3,  is  positive,  v.  > V.  > 0,  V.  > 0.75  ; 

3 3 - 3 3 ” 

(b)  Vc(1)  = " 7*56x10"3  ' v3<1)  = 4.23362. 

Clearly  u~(x;y  ) is  still  stable,  in  fact  (d/dx)  v7  ( 1;  y ) 

= 4.23  is  very  close  to  (d/dx) v^ (1; 1. 0)  = 4.31,  showing 

the  extremely  smooth  dependence  of  u^fxjy)  on  y. 

The  negativity  of  (d/dx)vc(l)  would  seem  to  point  Lo 

outright  instability  of  uc(x),  but  in  fact  its  small 

numerical  value  indicates  that,  within  the  accuracy  of  the 

computations,  it  is  actually  zero,  i.e.,  u (x)  is  neutrally 

c 

stable.  Indeed,  the  mathematical  theory  of  nonlinear 
problems  (Nirenberg,  1974)  shows  that  for  a bifurcation  to 
exist,  as  in  the  case  at  hand,  the  linearization 


LC  = • an  N (u,uc> Iu=u 


of  the  spatial,  time-independent,  operator  N at  the  bifurca- 
tion point  (u  , y ) has  to  have  a zero  eigenvalue.  This 
c c 

follows  from  the  infinite-dimensional  generalization  of  the 
one -dimensional  fact  that  a function  can  be  inverted  in  a 
neighborhood  of  a point  at  which  it  has  a nonzero  derivative , i.e., 
that  it  has  a unique  branch  near  such  a point. 

We  also  computed  the  lowest  eigenvalues  an<^ 

corresponding  eigenfunction  of  L (y  ) for  j = c,3, 

J J C 
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by  the  shooting  method  described  in  Subsection  5b.  In  this 
computation,  the  results  on  v ^ mentioned  before  were  valuable 
in  making  a first  guess  for  the  shooting  parameter  A.  The 
computations  yielded 


A(1)  = - 1.287xl0_11  , = 5.07265*1(T9  . 

c 3 

(i  j -9 

Again  we  notice  that  A^  (Uc)  = 5.07x10  is  very  close  to 
A^1)  (1.0)  = 5.13X10-9,  whereas  | A<5.1>  | = 0.4xl(f2  xj1)(1.0) 
= 0.2xio  A^  (1.0)  is  practically  zero,  as  it  has  to  be 


analytically . 

It  is  clear  by  the  continuity  of  (U)  in  the  parameter 

y that  for  y < y <.  1.05  we  have  x|^  > ^3^  > an<^ 

A^1^  <0,  so  that  the  interglacial  and  the  ice-covered  earth 
are  stable  for  the  entire  range  of  y explored,  whereas  the 
glacial  is  unstable  for  the  same  range  of  y.  Furthermore  the 
ice-covered  earth  is  stable  also  for  smaller  y. 

There  is  one  further  point  of  view,  which,  while  illumi- 
nating the  significance  of  the  neutral  stability  of  uc(x), 
also  argues  for  our  linear  stability  analysis  being  sufficient 
for  concluding  on  nonlinear  stability  or  instability  of  the 
steady-state  solutions  of  (4')  corresponding  to  different 
values  of  y.  This  viewpoint  has  to  do  with  the  existence  of 
a variational  principle  for  (4').  Indeed 

N (u; y)  = 0 

is  the  Euler-Lagrange  equation  for  the  extrema  of  the  functional 


J(u;y)  = 


pux  + rG(x,u)  f dx  ; 


here  p,  r are  given  by  (17b,c)  and 
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G(x,u) 


F(x,w)  doj  , 

where  F(x,u)  is  defined  by  (7d). 

Clearly  the  stable  solutions  u^(x;l),  u3(x;l)  corres- 
pond to  local  minima  of  J(u;l),  whereas  u2(x;l)  is  a local 
maximum.  As  u^(x;y),  which  is  a minimum  for  y > yc  , 
coallesces  with  u2(x;y),  which  is  a maximum  for  y > yc  , at 
y = yc  , a saddle  point  u = uc(x)  obtains,  whereas  u3(x;yc) 
is  still  a minimum. 

This  variational  interpretation  makes  it  very  plausible 
that  solutions  u(x,t;y)  of  the  "flow" 

ufc  = N (u; y)  , y > yc  , 

with  initial  conditions  near  u^ (x; y) , that  is  at  a finite 

but  small,  rather  than  infinitesimal,  distance  from  u^(x;y), 

j = 1,2,3,  will  tend  as  t -*■  00  towards  u^  if  u..  is  a minimum, 

i.e.,  j = 1,3,  and  away  from  it  when  u^  is  a maximum,  i.e., 

j = 2.  Similarly,  for  y = yc  , solutions  starting  near 

u,(x;y  ) will  still  converge  to  u-,  , but  solutions  starting 
3 c J 

near  uc(x) , though  they  may  hover  for  a long  time  near  uc  , 
since  x . , *►  °°  as  y y , will  eventually  move  away  from  it, 

along  a negative  slope  of  the  saddle,  and  finally  tend  towards 
the  absolute  minimum  u^.  This  seems  a rather  satisfactory, 
although  heuristic,  explanation  of  the  results  of  the  time- 
dependent  computations  of  Dwyer  and  Petersen  (1973)  and  of 
Schneider  and  Gal-Chen  (1973),  which  we  mentioned  alreay  in 
Section  5. 
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7.  Concluding  Remarks 


We  studied  the  zonally-and-vertically  averaged  energy- 
balance  climate  model  governed  by  equations  (l-3)  ; these 
equations  are  based  on  simple  parameterizations  of  albedo, 
greenhouse  effect  and  eddy  diffusion  of  heat  in  terms  of 
yearly  averaged  sea-level  temperature,  which  is  the  only 
dependent  variable  of  the  model. 

Three  positive  steady- state  solutions  of  the  model, 
symmetric  with  respect  to  the  equator,  were  found  by 
accurate  numerical  computations , and  apparently  no  more 
such  solutions  exist.  These  steady  states  can  be  identified 
with  an  interglacial  climate,  approximating  very  well  the 
one  prevailing  presently  on  earth,  a glacial  climate,  and 
a climate  during  which  the  earth  would  be  completely  ice 
covered.  The  climates  obtained  were  only  slightly  changed 
when  making  small  changes  in  the  numerical  values  of  the 
coefficients  and  when  making  certain  changes  in  the  functional 
form  of  the  model's  equations.  However,  the  bounds  on  the 
values  the  albedo  can  take  were  essential  in  order  to  obtain 
these  three  climates;  also  linearizing  the  outgoing  planetary 
radiation  resulted  in  a reduction  of  the  number  of  solutions. 

We  then  determined  the  stability  of  the  time-dependent 
solutions  of  (1')  under  small  perturbations  about  the  model's 
steady  states.  This  stability  was  shown  to  depend  on  certain 
properties  of  a comparison  function,  which  was  constructed 


numerically, 
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We  found  that  the  interglacial  and  the  "deep  freeze" 
climate  are  stable,  and  that  the  glacial  climate  is  unstable. 
This  means  that  the  first  two  can  obtain,  at  least  approximate- 
ly, as  steady  states  in  a physical  system  governed  by  equa- 
tions very  similar  to  (1-3),  but  that  the  latter  cannot; 
the  same  is  true  about  these  climates  as  limiting  steady 
states  for  time -dependent  numerical  solutions  of  such 
equations . 

We  further  showed  how  changes  in  an  important  parameter, 
the  average  intensity  of  the  solar  radiation,  influence  the 
steady-state  solutions  of  the  model.  The  dependence  on  this 
parameter  of  all  steady  states  was  shown  to  be  gradual  and 
smooth  for  increases  of  up  to  5%  and  decreases  of  up  to 
about  2%.  However  for  a critical  value  of  the  parameter, 
equal  to  98.15%  of  its  present  value,  the  glacial  and  inter- 
glacial climates  coallesced  and  they  disappeared  entirely  for 
smaller  values  of  the  parameter,  leaving  the  ice-covered  earth 
as  the  only  possible  stable,  steady  climate  of  the  model. 

This  result  is  important,  as  it  stresses  the  difference 
between  the  stability  of  a steady  state  with  respect  to  the 
time  evolution  of  a physical  system  governed  by  a given, 
fixed  equation,  and  the  stability  of  a steady  state  with 
respect  to  changes  in  a parameter,  which  determines  the 
behavior  of  the  system.  For  definiteness,  let  us  call  the 
former  internal  stability  and  the  latter  external  stability; 
we  have  shown  that  the  "deep  freeze"  is  stable  for  our  system 

-48- 


MM 


■r Ml-OIIH  >.  * "BfanHOM  *P  *•  * 


both  internally  and  externally,  that  the  glacial  is  unstable 
in  both  senses,  and  that  the  interglacial,  or  "present 
climate,"  is  internally  stable,  but  externally  unstable. 

The  limitations  of  equations  (1-3)  as  a model  for  the 
description  of  the  long-term  behavior  of  the  atmosphere- 
ocean-cryosphere  system,  and  of  energy-balance  models 
in  general,  have  been  discussed  extensively.  Because  of 
these  limitations,  we  believe  that  the  results  above 
should  not  be  taken  at  face  value  as  statements  about  the 
climate  of  our  planet.  These  results,  however,  seem  to 
clarify  the  physical  content  and  mathematical  properties 
of  such  models.  Also,  the  methods  used  here  could  be 
helpful  in  investigating  other  models,  which  will  include 
more  elaborate  and  reliable  parameterizations  of  the 
physical  phenomena  governing  climate.  We  further  hope  that 
insight  gained  into  the  behavior  of  solutions  of  a certain 
type  of  model  will  advance  the  formulation  of  other  models, 
and  that  these  will  come  closer  to  explaining  past  changes 
in  climate  and  predicting  future  changes. 
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Table  Captions 


Table  1.  Empirical  functions  appearing  in  equation  (4) . 

The  functions  Q,  b,  z are  based  on  data  in  Tables  1 
and  2 of  Sellers  (1973).  The  functions  C,  , k 2 
are  based  on  data  provided  by  Dr.  T.  Gal-Chen  (1974, 
personal  communication) , and  used  in  the  (SV)  model 
of  Schneider  and  Gal-Chen  (1973). 

Table  2.  Influence  of  different  modifications  in  the  model's 
equation  (4)  on  the  number  of  steady-state  solutions 
and  the  numerical  values  of  these  solutions.  The 
existing  solutions  are  identified  by  the  temperature 
at  the  pole,  u^  (0)  , j = 1,2,3.  In  case  a solution  is 
missing,  this  is  indicated  by  (-)  in  the  corresponding 
row- and- column  location.  S stands  for  the  coefficients 
being  fitted  by  cubic  splines,  B for  Bernstein  poly- 
nomials. A downward  arrow  (+)  to  the  right  of  a comment 
indicates  that  the  equation  used  in  the  numerical 
experiments  reported  in  all  subsequent  rows  had  the 
feature  pointed  out  in  that  comment.  Otherwise  comments 
refer  only  to  the  row  in  which  they  appear.  A left 
arrow,  x «-  y,  means  that  the  quantity  x was  replaced 
in  the  equation  by  the  quantity  y.  The  entries  given 
with  less  than  five  decimal  digits  resulted  from  computa 
tions  with  lower  precision  than  indicated  in  the  text. 
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Table  3.  Dependence  of  the  steady  states  (x;  p)  ,1^  (x;  y)  on  the 
normalized  average  intensity  of  the  solar  radiation,  u. 

The  columns  give  u(0),  the  temperature  at  the  pole, 

6^u,  the  average  temperature,  and  Au  = u(l)  - u(0), 
the  pole- to-equator  temperature  difference  for  u^  and 
for  U2  respectively.  The  values  for  p = 0.98151822 
correspond  to  the  bifurcating  steady  state  u = u (x) . 
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Figure  Captions 


Figure  1 

Comparison  of  curve  fitting  by  (i)  Bernstein  poly- 
nomial approximation,  indicated  by  a dash-dot  line,  and  by 
(ii)  cubic  spline  interpolation,  indicated  by  a solid  line. 
Bernstein  polynomials  are  not  in terpolatory  and  they  are 
variation  diminishing,  i.e.,  they  have  the  property  of 
smoothing  out  the  data;  this  results  in  a rather  poor 
approximation.  Cubic  splines  are  not  variation  diminishing 
and  they  are  very  good  approximants . 

Figure  la 

For  a very  smooth  function,  like  u(x),  the  two  approxi- 
mation procedures  yield  curves  very  close  to  each  other. 
Figure  lb 

For  a function  of  large  total  variation,  like  k(x,u(x)), 
the  two  procedures  yield  curves  which  can  differ  pointwise  by 
as  much  as  50  percent  of  the  average  value. 

Fi gure  2 

Numerically  obtained  values  of  ux(l;Ug)  as  a function 

Of  Ug  = u (0) . 

Figure  2a 

Comparison  of  the  results  for  equations  (7),  in 
which  k = k(x,u),  with  those  for  equations  (!'),  in  which 
k = k(x,u(x));  the  results  for  (7)  are  indicated  by  a solid 
line,  those  for  (7')  by  a dash-dot  line. 
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Figure  2b 

Results  of  (7')  for-  -1330K  < u(0)  £ 300  K.  Notice  that 
as  Uq  = u(0)  tends  towards  the  ends  of  the  interval, 
ux(l;uQ)  + +°°.  The  solution  u4(x),  corresponding  to  the 
negative  root  of  this  curve,  Uq  - -186  K,  does  not  have  a 
physical  significance. 

Figure  3 

Comparison  of  the  solutions  of  (4),  indicated  by  a solid 
line,  with  those  of  (4')»  indicated  by  a dash-dot  line.  The 
circles  indicate  mesh  data  for  u = u(x) . 

Figure  3a 

Values  of  the  solutions  u^ (x) , j = 1,2,3,  for  (4)  and 
for  (4').  The  respective  values  for  j = 1,3  are  practically 
indistinguishable,  whereas  for  j = 2 a slight  difference 
exists  between  the  solution  of  (4)  and  that  of  (4  ) . 

Figure  3b 

Values  of  the  derivatives  ^ u^(x),  j = 1,2,3,  for  (4) 
and  for  (4').  The  differences  are  larger  than  in  the  function 

values  themselves. 


Figure  4 


.(1) 


The  first  eigenfunctions,  wj"'  ^ • 3 - 1*2,3  , of 


Lj  - ' fe  N(u)'u=u. 


Figure  5 


The  function  u (1;Uq),  obtained  by  integrating  (7  ) 


numerically  with  y = yc  and  with  different  values  of  Uq, 
u~  = u(0)  > 120  K. 


Figure  6 

Dependence  of  the  solutions  Uj(x;p),  j = 1,2,  on  the 
parameter  p.  The  two  plots  for  Uj(0),  are  very  simi- 

lar; the  plot  for  Au^  = Uj(l)-Uj(0)  is  rather  different, 
although  it  exhibits  the  same  behavior  in  the  neighborhood 
of  the  critical  point  c.  The  circles  indicate  the  values 
actually  computed,  for  p = pc  , 0.982,  0.985,  0.99,  1.00, 

1.01,  1.02,  1.03,  1.04,  1.05.  The  letter  c distinguishes 
the  values  of  the  plotted  quantities  u ^ ( 0 ) , £^u^  , Au^ 

corresponding  to  the  bifurcating  solution  u = uc(x). 

Figure  7 

The  bifurcating  solution  u = uc(x).  Notice  that  the 
ice  line,  which  corresponds  to  u = 273  K,  is  at  about  45°  lat. 
for  this  solution,  i.e.,  an  ice  cover  extending  beyond  this 
latitude  would  eventually  cover  the  entire  earth. 
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120.  140.  160.  180.  200.  220.  240.  260.  280.  300.  U(0) 

fig.  5 
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